- Structure of spatial data in R
- Reading / Writing Spatial Data
- Visualizing Spatial Data
- Analysing Spatial Data
Wednesday, February 18 2015
| Mandatory Components | Optional Components |
|---|---|
| .shp - actual geometry of feature | .prj - CRS and projection info in WKT format |
| .shx - shape index - binary file giving position index | .sbn and .sbx- spatial indexing files |
| .dbf - attribute information, stored in dBase IV format | .xml - metadata file |
library(rgdal) data(nor2k) class(nor2k)
[1] "SpatialPointsDataFrame" attr(,"package") [1] "sp"
slotNames(nor2k)
[1] "data" "coords.nrs" "coords" "bbox" "proj4string"
plot(nor2k, axes=T,
main='Peaks in Norway over 2000 meters')
library(rgdal) data(nor2k) summary(nor2k)
Object of class SpatialPointsDataFrame
Coordinates:
min max
East 404700 547250
North 6804200 6910050
Height 2001 2469
Is projected: TRUE
proj4string :
[+proj=utm +zone=32 +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0]
Number of points: 300
Data attributes:
Nr. Navn Kommune
Min. : 1.00 Length:300 Length:300
1st Qu.: 75.75 Class :character Class :character
Median :150.50 Mode :character Mode :character
Mean :150.50
3rd Qu.:225.25
Max. :300.00
library(raster) r <- raster(ncol=10, nrow=10) values(r) <- c(1:100) plot(r, main='Raster with 100 cells')
library(maps);library(sp);require(knitr) data(us.cities) knitr::kable(us.cities[1:5,])
| name | country.etc | pop | lat | long | capital |
|---|---|---|---|---|---|
| Abilene TX | TX | 113888 | 32.45 | -99.74 | 0 |
| Akron OH | OH | 206634 | 41.08 | -81.52 | 0 |
| Alameda CA | CA | 70069 | 37.77 | -122.26 | 0 |
| Albany GA | GA | 75510 | 31.58 | -84.18 | 0 |
| Albany NY | NY | 93576 | 42.67 | -73.80 | 2 |
class(us.cities) # simple data frame
[1] "data.frame"
library(maps);library(sp)
data(us.cities)
coordinates(us.cities) <- c("long", "lat")
class(us.cities)
[1] "SpatialPointsDataFrame" attr(,"package") [1] "sp"
plot(us.cities, pch = 20, col = 'forestgreen', axes=T,
xlim=c(-125,-70), ylim=c(26,55))
library(maps) par(mfrow=c(1,1)) map()
map.text('county','oregon')
map.axes()
title(main="Oregon State")
Loading administrative backgrounds from Global Administrative Areas is another good option (http://gadm.org/)
library(rgdal)
setwd('K:/GitProjects/RUserWebinar')
HUCs <- readOGR('.','OR_HUC08')
writeOGR(HUCs, '.','HUC', driver='ESRI Shapefile')
library(raster)
r <- raster('clay.tif')
# crop it
r <- crop(r, extent(-1000000, 2000000, 100000, 3000000) )
writeRaster(r, 'clay_smaller.tif',format='GTiff')
require(sp);require(rgeos);load("K:/GitProjects/RUserWebinar/Data.RData")
# A spatial PolygonsDataFrame has 5 slots:
str(HUCs, 2)
Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots ..@ data :'data.frame': 91 obs. of 3 variables: ..@ polygons :List of 91 ..@ plotOrder : int [1:91] 78 1 32 87 81 24 26 11 38 3 ... ..@ bbox : num [1:2, 1:2] -2297797 2191569 -1576981 2914667 .. ..- attr(*, "dimnames")=List of 2 ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
require(sp);require(rgeos);load("K:/GitProjects/RUserWebinar/Data.RData")
# Each polygon element has 5 of it's own slots - here we look at first one:
str(HUCs@polygons[[1]])[]
Formal class 'Polygons' [package "sp"] with 5 slots ..@ Polygons :List of 1 .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots .. .. .. ..@ labpt : num [1:2] -1805054 2272669 .. .. .. ..@ area : num 9.13e+09 .. .. .. ..@ hole : logi FALSE .. .. .. ..@ ringDir: int 1 .. .. .. ..@ coords : num [1:13339, 1:2] -1782485 -1782433 -1782329 -1782239 -1782084 ... ..@ plotOrder: int 1 ..@ labpt : num [1:2] -1805054 2272669 ..@ ID : chr "0" ..@ area : num 9.13e+09
NULL
require(sp);require(rgeos); load("K:/GitProjects/RUserWebinar/Data.RData")
# Here we get a vector of area for features in HUCs spdf
area(HUCs[1:4,])
[1] 9128424113 2997942053 6065238017 4861449846
# Total Area - gArea function in rgeos gives same result sum(area(HUCs))
[1] 312620322426
# Area of a particular feature HUCs@polygons[[1]]@Polygons[[1]]@area
[1] 9128424113
require(sp);load("K:/GitProjects/RUserWebinar/Data.RData")
plot(HUCs, axes=T, main='HUCs in Oregon')
require(sp)
load("K:/GitProjects/RUserWebinar/Data.RData")
plot(HUCs, axes=T, main='HUCs in Oregon')
# Function to calculate percent of area
AreaPercent <- function(x) {
tot_area <- sum(sapply(slot(x, "polygons"),
slot, "area"))
sapply(slot(x, "polygons"), slot,
"area") / tot_area * 100
}
# just plot bigger HUCs
plot(HUCs[AreaPercent(HUCs) > 1,], add=T,
col='Steel Blue')
require(ggplot2)
load("K:/GitProjects/RUserWebinar/Data.RData")
# Take a look at some USGS stream gages for PNW
gages@data[1:5,5:8]
## LON_SITE LAT_SITE NearCity FLOW ## 1 -123.3178 42.43040 Grants Pass 3402.544 ## 2 -122.6011 42.42763 Eagle Point 60.201 ## 3 -119.9233 42.42488 Altamont 34.272 ## 4 -122.6011 42.40819 Eagle Point 104.177 ## 5 -122.5373 42.40263 Eagle Point 72.511
# Explicitly set CRS for layers
gages <- spTransform(gages,
CRS('+init=epsg:2991'))
# Locations of gages
ggplot(gages@data, aes(LON_SITE, LAT_SITE)) +
geom_point(aes(color=log10(FLOW))) +
coord_equal()
load("K:/GitProjects/RUserWebinar/Data.RData")
gages_proj <- spTransform(gages,
CRS('+init=epsg:2991'))
HUCs_proj <- spTransform(HUCs,
CRS('+init=epsg:2991'))
HUCs_proj$LON <- coordinates(HUCs_proj)[,1]
HUCs_west <- HUCs_proj[HUCs_proj$LON < 400000, ]
options(scipen=3)
plot(gages_proj, axes=T, col='blue')
plot(HUCs_west, add=T)
load("K:/GitProjects/RUserWebinar/Data.RData")
gages_west <- gages_proj[HUCs_west,]
plot(HUCs_west, axes=T)
plot(gages_west, add=T, col='blue')
load("K:/GitProjects/RUserWebinar/Data.RData")
OverUpdate <- function(points, polys) {
pointpoly <- over(points, polys)
points@data <- data.frame(points@data, pointpoly)
}
gages_proj@data <- OverUpdate(gages_proj, HUCs_proj)
head(gages_proj@data[,c(3,12)])
STATION_NM HUC_8 1 ROGUE RIVER AT GRANTS PASS, OR 17100308 2 N F LTL BUTE CR AB INTKE CANL LKECREEK OREG 17100307 3 HONEY CREEK NEAR PLUSH,OREG. 17120007 4 SOUTH FORK LITTLE BUTTE CR NR LAKECREEK,OREG. 17100307 5 NO FK LITTLE BUTTE CR NR LAKECREEK,OREG. 17100307 6 TWELVEMILE CREEK NEAR PLUSH,OREG. 17120007
load("K:/GitProjects/RUserWebinar/Data.RData")
HUCs_proj$StreamFlow <- over(HUCs_proj,gages_proj[8],fn = mean)
head(HUCs_proj@data[!is.na(HUCs_proj@data$StreamFlow),])
HUC_8 SUM_ACRES LON FLOW 2 17050103 1498689 716737.1 32.78000 5 17050107 956848 669298.4 930.73000 8 17050110 1264239 631925.6 622.90075 10 17050116 1554027 572126.3 130.40100 11 17050117 607018 633179.4 304.05640 12 17050118 375002 617941.8 31.72333
load("K:/GitProjects/RUserWebinar/Data.RData")
head(cities[,c(3:4, 7:8)])
NAME ST_ABBREV POP_2000 POP2007 1 Astoria OR 9813 9901 2 Warrenton OR 4096 4413 3 Gearhart OR 995 1033 4 Seaside OR 5900 5982 5 Clatskanie OR 1528 1664 6 Cannon Beach OR 1588 1669
require(plyr);load("K:/GitProjects/RUserWebinar/Data.RData")
# We can use match or join to connect to our spatial gages
# gages$POP2007 <- cities$POP2007[match(gages$NearCity, cities$NAME)]
# OR set a common name field and use join from plyr
names(gages)[7] <- "NAME"
gages@data <- join(gages@data, cities)
head(gages@data[,c(3,18)])
STATION_NM POP2007 1 ROGUE RIVER AT GRANTS PASS, OR 24753 2 N F LTL BUTE CR AB INTKE CANL LKECREEK OREG 6819 3 HONEY CREEK NEAR PLUSH,OREG. 20493 4 SOUTH FORK LITTLE BUTTE CR NR LAKECREEK,OREG. 6819 5 NO FK LITTLE BUTTE CR NR LAKECREEK,OREG. 6819 6 TWELVEMILE CREEK NEAR PLUSH,OREG. 20493
library(rgeos); library(rgdal)
lps <- coordinates(HUCs)
IDFourBins <- cut(lps[,1], quantile(lps[,1]),
include.lowest=TRUE)
regions = gUnaryUnion(HUCs,IDFourBins)
regions = SpatialPolygonsDataFrame(regions,
data.frame(regions = c('Coastal',
'Mountains','High Desert','Eastern')),
match.ID = FALSE)
plot(regions, axes=T)
text(coordinates(regions),
label = regions$regions, cex=.8)
require(raster)
load("K:/GitProjects/RUserWebinar/Data.RData")
clay <- raster('clay.tif')
inMemory(clay)
[1] FALSE
cellStats(clay, min); cellStats(clay, max)
[1] 0
[1] 73.6665
projection(clay)
[1] "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
plot(clay)
require(raster)
clay <- raster('clay.tif')
clay_OR <- crop(clay, extent(-2261000, -1594944, 2348115, 2850963))
plot(clay_OR)
require(raster)
clay <- raster('clay.tif')
gages <- spTransform(gages, CRS(projection(clay)))
gages$clay = extract(clay,gages)
head(gages@data[,c(3,12)])
STATION_NM clay 1 ROGUE RIVER AT GRANTS PASS, OR 24.0840 2 N F LTL BUTE CR AB INTKE CANL LKECREEK OREG 37.8365 3 HONEY CREEK NEAR PLUSH,OREG. 28.0060 4 SOUTH FORK LITTLE BUTTE CR NR LAKECREEK,OREG. 37.8365 5 NO FK LITTLE BUTTE CR NR LAKECREEK,OREG. 37.8365 6 TWELVEMILE CREEK NEAR PLUSH,OREG. 40.8730
library(rasterVis)
alt <- getData('worldclim', var='alt', res=2.5)
usa1 <- getData('GADM', country='USA', level=1)
oregon <- usa1[usa1$NAME_1 == 'Oregon',]
alt <- crop(alt, extent(oregon) + 0.5)
alt <- mask(alt, oregon)
levelplot(alt, par.settings=GrTheme)
require(plotGoogleMaps)
HUCs <- spTransform(HUCs, CRS('+init=epsg:28992'))
map <- plotGoogleMaps(HUCs, filename='RGoogleMapsExample.htm')
library(ggmap) mymap <- get_map(location = "Corvallis, OR", source="google", maptype="terrain",zoom = 12) ggmap(mymap, extent="device")